不得不爱（1）——网状meta分析必备技能（三）：R软件gemtc包
我们仅仅是代码的编辑者、整合者、搬运工，仅免费传授方法，下文数据和代码取自于网络和免费软件“R语言说明书”，如果您觉得我们侵犯了您的版权，请通知我们撤稿。
我们没有详细讲解理论是因为——高手们已经发表了网状meta分析的理论文章，请大家在cnki检索“网状meta分析”学习理论知识（高手告诉我：复制超过6个字，意味着侵犯版权）。请大家谅解，谢谢！
###FS的宗旨：科学自由分享、人人平等，共求真理###
第一部分 二分类数据网状meta分析的异质性###########
3.1 新的一天了，大家一定很累了，不知道上一章（忘记了？点这里），大家操练好了没有呢？我们来点简单的baidu google……
今天在这里我们将“R软件gemtc包”的实战讲透彻！具体理论请自行学习，可以baidu gemtc，摸索一下。或者google科学上网。
或者使用猴哥为大家准备的gemtc官方的说明书，具体下载地址如下：
http://pan.baidu.com/s/1dFNUNyp
3.2 0001 二分类数据 残差 残差分布
3.2.1 残差 残差分布
在回归分析中，测定值与按回归方程预测的值之差，以δ表示。残差δ遵从正态分布N(0，σ2)。（δ-残差的均值）/残差的标准差，称为标准化残差，以δ*表示。δ*遵从标准正态分布N(0，1)。实验点的标准化残差落在(-2，2)区间以外的概率≤0.05。若某一实验点的标准化残差落在(-2，2)区间以外，可在95%置信度将其判为异常实验点，不参与回归直线拟合。
显然，有多少对数据，就有多少个残差。残差分析就是通过残差所提供的信息，分析出数据的可靠性、周期性或其它干扰。
3.2.2 代码如下，请复制到Rstudio的r软件中
猴哥为大家准备的苹果版和WINDOWS R软件和stata软件，具体下载地址如下：
http://pan.baidu.com/s/1dFNUNyp
#1清除环境变量，设置路径--------------------------------------------------
rm(list=ls())
getwd()
# 猴哥是放在F盘的0108gemtc_sm文件夹下，你可以自己新建一个文件夹
setwd("F:\\0108gemtc_sm\\gemtc\\inst")
# 如果没有安装gemtc，请执行如下命令，执行时去掉#，并在Rstudio中按 ctrl+enter
# install.packages("gemtc")
#2载入程序包--------------------------------------------------
# install.packages('code')
library(gemtc)
# library(code)
library(codetools)
library(testthat)
#3导入数据load file 1
#3.1第一种方式
导入执行包gemtc包中的系统文件，并命名为file
file <- system.file('extdata/luades-smoking.gemtc', package='gemtc')
# 看一下文件位置
print(file)
##读取gemtc文件，并保存为net
net = read.mtc.network(file)
# 或者命名为network
network <- read.mtc.network(file)
# 看一下文件数据
print(net)
print(network)
#3.2保存数据
# 这里net 和network是list，所以不能直接保存数据，我们只能保存数据框
write.csv(network$treatments, file = "F:\\0108gemtc_sm\\network0001_treatments.csv")
write.csv(network$data.ab, file = "F:\\0108gemtc_sm\\network0001_data.ab.csv")
# #或者，这里的file = " F:\\0108gemtc_sm\\network0001.csv"="F:/0108gemtc_sm/net0001.csv"
write.csv(net$treatments, file = "F:/0108gemtc_sm/net0001_treatments.csv")
write.csv(net$data.ab, file = "F:/0108gemtc_sm/net0001_data.ab.csv")
#这里说明了戒烟的4个治疗手段，A B C D
#这里说明了戒烟的24个研究，共计50个arm的情况
#4导入数据,第2种方法----------------------------------------
#大家习惯了在excel中处理数据，可以先转化为 .csv格式
# 读取CSV文件
# header=T,表示文件存在，表示需要读取数据的数据头部
# skip=0,表示数据中没有需要跳过的行, ./treatments0001.csv表示本文件夹中的treatments0001.csv文件，row.names = 1是第一列作为行名
treatments<- read.csv("./net0001_treatments.csv",header=T,sep = ",",skip=0,row.names = 1)
# 我们看一下数据
head(treatments)
# 我们看一下数据
data <- read.csv("./net0001_data0001.csv",header=T,sep = ",",row.names = 1)
head(data)
#5创建network对象，建立成network 单臂长数据格式,
description=" Luades_smoking "是对network的描述
network<- mtc.network(data, description="Luades_smoking", treatments=treatments)
#network是单臂长数据
network
# 进一步了解 network的属性
head(network)
#6构建模型####
#一致性模块
# 其中， network 为 network 数据， type 为是否选取一致性模型，
# n.chain 为迭代运算中链的条数。需要注意的是，目前该程序包通过 R 软件只能使用一致性模型
# model <- mtc.model(network, linearModel='fixed', likelihood='normal', link='identity')
# model1
model <- mtc.model(network, type = 'consistency')
# 以下model2模型和上面的模型，意义差不多
# model2
model <- mtc.model(network,type = "consistency", factor = 2.5, n.chain = 3)
#7网状证据图####
#网状图在可这里可以直接出图
#模型网状图
plot(network)
# # points(network, cex = .5, col = "dark red") #改变颜色没有用
# #network图汇总
summary(network)
# 网状图导出
#绘制tiff图片并保存在Rwork文件夹下
tiff(file="network.tiff")
plot(network)
dev.off()
#8运行结果####
#使用#JAGS进行运算
#8.1简单命令
result <- mtc.run(model)
#8.2复杂命令
result <-mtc.run(model, sampler ="rjags", n.adapt = 5000, n.iter = 20000, thin = 1)
#Setting the sampler is deprecated, only JAGS is supported. #sunmin dell PC
#9拟合模型，等待一段时间
#####################以下内容和上一章相同，是结果的输出部分######################
#10网状证据图####
#网状图在可这里可以直接出图
#模型网状图
plot(network)
# # points(network, cex = .5, col = "dark red") #改变颜色没有用
# #network图汇总，右下角的窗口中出现了，我们要的图形
# 网状图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0101network.tiff")
plot(network)
dev.off()
#文件夹中出现了，我们要的森林图，其他同理可得
summary(network)
#诊断收敛性
gelman.plot(result)
# 收敛图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0102gelman.plot.tiff")
gelman.plot(result,auto.layout =F)
dev.off()
#结果汇总
summary(result)
#密度图
plot(result)
# 密度图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0103密度图.tiff")
plot(result)
dev.off()
#森林图
forest(result)
# 森林图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0104森林图.tiff")
forest(result)
dev.off()
#等级排名
ranks <- rank.probability(result)
print(ranks)
#排序图，堆积排序图
tiff(file="0105堆积排序图.tiff")
plot(ranks) # plot a cumulative rank plot
dev.off()
#单个排序图，等级图
tiff(file="0106单个排序图.tiff")
barplot(t(ranks), beside=TRUE) # plot a 'rankogram'
dev.off()
#单个排序图，等级图，我们让色彩丰富一些
tiff(file="0107单个排序图.tiff")
barplot(t(ranks), beside = TRUE,
col = c("lightblue", "mistyrose", "lightcyan",
"lavender", "cornsilk"),
legend = rownames(VADeaths), ylim = c(0, 1))
title(main = "Death Rates in Virginia", font.main = 8)
dev.off()
#单个排序图，等级图，我们让色彩更丰富一些
barplot(t(ranks), beside=TRUE, col = colors(1:14))
# 单个排序图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0108单个排序图.tiff")
barplot(t(ranks), beside=TRUE, col = colors(1:14))
dev.off()
#11相对影响森林图
# 相对影响森林图导出vsB
tiff(file="0109相对影响森林图导出vsA.tiff")
forest(relative.effect(result, "A"))
dev.off()
summary(relative.effect(result, "A", c("B", "C", "D")))
#12残差
print(result$deviance)
#13杠杆效应####
## residual deviance plot
# #the first method，暂不能实现，但需要运行
if (model$data$ns.r2 + model$data$ns.rm == 0) {
tpl <- gemtc:::arm.index.matrix(model[['network']])
study <- matrix(rep(1:nrow(tpl), times=ncol(tpl)), nrow=nrow(tpl), ncol=ncol(tpl))
study <- t(study)[t(!is.na(tpl))]
devbar <- t(results$deviance$dev.ab)[t(!is.na(tpl))]
title <- "Per-arm residual deviance"
xlab <- "Arm"
} else {
nd <- model$data$na
nd[-(1:model$data$ns.a)] <- nd[-(1:model$data$ns.a)] - 1
devbar <- c(apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE), result$deviance$dev.re) / nd
study <- 1:length(devbar)
title <- "Per-study mean per-datapoint residual deviance"
xlab <- "Study"
}
plot(devbar, ylim=c(0,max(devbar, na.rm=TRUE)),
ylab="Residual deviance", xlab=xlab,
main=title, pch=c(1, 22)[(study%%2)+1])
# #the second method
###########################
for (i in 1:length(devbar)) {
lines(c(i, i), c(0, devbar[i]))
}
# w <- sign(r - rfit) * sqrt(devbar)
# plot(w, leverage, xlim=c(-3,3), ylim=c(0, 4.5))
# x <- seq(from=-3, to=3, by=0.05)
# for (c in 1:4) { lines(x, c - x^2) }
fit.ab <- apply(result$deviance$fit.ab, 1, sum, na.rm=TRUE)
dev.ab <- apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE)
lev.ab <- dev.ab - fit.ab
fit.re <- result$deviance$fit.re
dev.re <- result$deviance$dev.re
lev.re <- dev.re - fit.re
nd <- model$data$na
nd[-(1:model$data$ns.a)] <- nd[-(1:model$data$ns.a)] - 1
w <- sqrt(c(dev.ab, dev.re) / nd)
lev <- c(lev.ab, lev.re) / nd
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
# 重新运行下画图程序，可以出图
tiff(file="0110相对影响森林图导出vsB.tiff")
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
dev.off()
#14万众期待的异质性图，猴哥自家的代码，google里面没有哦############
# mtc.anohe(network)
# 14非常重要##############
# 14smok研究间的异质性和森林图##############
result.anohe <- mtc.anohe(network, n.adapt=1000, n.iter=5000)
# plot(result.anohe)
summary.anohe <- summary(result.anohe)
plot(summary.anohe, xlim=log(c(0.2, 5)))
summary.anohe
# Analysis of heterogeneity
# =========================
#
# Per-comparison I-squared:
# -------------------------
#
# t1 t2 i2.pair i2.cons incons.p
# 1 A B 61.68561 66.88978 0.7106796
# 2 A C 92.68188 92.19784 0.9117685
# 3 A D 86.36661 57.96352 0.4816269
# 4 B C 0.00000 0.00000 0.4877663
# 5 B D 40.96877 0.00000 0.8934959
# 6 C D 62.36553 58.46978 0.2931933
#
# Global I-squared:
# -------------------------
#
# i2.pair i2.cons
# 1 89.69921 88.5369
summary.anohe$consEffects
# t1 t2 pe ci.l ci.u
# 1 A B 0.4815935 -0.2926061 1.315680
# 2 A C 0.8316086 0.3862057 1.361283
# 3 A D 1.0858208 0.2467852 2.027102
# 4 B C 0.3514423 -0.4599256 1.200236
# 5 B D 0.6080877 -0.3433846 1.603331
# 6 C D 0.2507416 -0.5646447 1.088071
summary.anohe$studyEffects
# study t1 t2 pe ci.l ci.u
# 1 1 A C 1.08468247 0.131825883 2.03753906
# 2 1 A D 0.13179134 -1.246876317 1.51045900
# 3 1 C D -0.95289113 -1.864452746 -0.04132951
# 4 2 B C 0.01477717 -1.344989257 1.37454359
# 5 2 B D 0.26425258 -0.598264269 1.12676944
# 6 2 C D 0.24947542 -0.590747993 1.08969883
# 7 3 A B -0.01509682 -0.346692546 0.31649890
# 8 4 A B 0.39503472 -0.244562172 1.03463162
# 9 5 A B 0.72634581 -0.155105472 1.60779709
# 10 6 A C 2.21424725 1.930107850 2.49838665
# 11 7 A C 1.09410003 -0.675247258 2.86344732
# 12 8 A C 0.42748862 0.122258254 0.73271898
# 13 9 A C 18.71338529 -8.494560712 45.92133130
# 14 10 A C 2.86935813 1.584637810 4.15407844
# 15 11 A C 3.08048940 0.444298661 5.71668015
# 16 12 A C 0.50230300 -0.532874788 1.53748079
# 17 13 A C 0.46393989 0.179475792 0.74840399
# 18 14 A C -0.13551085 -0.764993165 0.49397147
# 19 15 A C -0.24361822 -0.585177256 0.09794082
# 20 16 A C 0.04031725 -0.329491827 0.41012633
# 21 17 A C 0.39445062 0.062293725 0.72660752
# 22 18 A C 0.14917935 -1.067970279 1.36632897
# 23 19 A C 0.58502492 -0.002707489 1.17275734
# 24 20 A D 24.42541723 -2.919989916 51.77082438
# 25 21 B C -0.15718293 -1.008726421 0.69436056
# 26 22 B D 1.10500752 0.174982175 2.03503287
# 27 23 C D 0.70010819 -0.112046095 1.51226247
# 28 24 C D -0.51648572 -1.993457795 0.96048635
summary.anohe$pairEffects
# t1 t2 pe ci.l ci.u
# 1 A B 0.34239499 -0.71419051 1.465357
# 2 A C 0.82173062 0.33285861 1.368663
# 3 A D 1.65193949 -0.04405122 3.577416
# 4 B C -0.08185359 -1.55575229 1.402015
# 5 B D 0.67841016 -0.72818479 2.105401
# 6 C D -0.07257150 -1.12055323 0.925233
network
result <-mtc.nodesplit(network)
summary(result)
names(result)
# [1] "d.A.C" "d.A.D" "d.B.D" "d.C.D" "consistency"
summary(result$d.A.C)
# Overall summary and plot
summary.ns <- summary(result)
print(summary.ns)
plot(summary.ns)
# 建议大家在 r软件中画图哦，Rstudio中只能整理看一下结果
#或者直接导出到pdf
###FS的宗旨：科学自由分享、人人平等，共求真理###
中文题目：2016rNMA gemtc五种常见麻醉在小儿麻醉和恢复中的特点：网络的Meta分析
Emergence and Recovery Characteristics of Five Common
Anesthetics in Pediatric Anesthesia: a Network Meta-analysis
Jianrong Guo1 & Xiaoju Jin2 & Huan Wang1 & Jun Yu2 & Xiaofang Zhou1 & Yong Cheng1 &Qiang Tao1 & Li Liu1 & Jianping Zhang Mol Neurobiol
DOI 10.1007/s12035-016-9982-3 （在对话框回复NMA01，获得文献哦）
图1 是r软件gemtc包做的5种治疗方案的 网状证据图 代码见我们上面的 “#7网状证据图####”部分
图2 相对相应森林图
代码见我们的“# 森林图导出” 和 “#11相对影响森林图”
图3和图5 是探讨模型的一致性和不一致性，代码见我们上文的“#15节点劈裂法,探讨模型的一致性和不一致性####“
图4和图6 是探讨模型的一致性和不一致性，代码见我们上文的“#15节点劈裂法,探讨模型的一致性和不一致性####“
图7 是探讨模型的一致性和不一致性，代码见我们上文的“#15节点劈裂法,探讨模型的一致性和不一致性###”
图8 是探讨模型的一致性和不一致性，代码见我们上文的“#15节点劈裂法,探讨模型的一致性和不一致性###”
###FS的宗旨：科学自由分享、人人平等，共求真理###
谢谢大家，请多指教！！！
猴哥：Freescience公众号meta分析栏目现任主编。副主任医师，武汉大学肿瘤学博士，专注于胃肠道肿瘤分子生物学机制、系统评价/Meta分析、数据挖掘、临床统计研究。
科研路，不孤单！
Freescience医学科研联盟全国火热招募ing
50家高校及医院的小伙伴已经加入啦
具体点这里
想围观一篇SCI论文是怎样写成的吗？
Freescience线上沙龙之SCI论文合作呼唤你
具体点这里
FS科研软件库，集合50+医学科研必备神器，现在统统打包分享，点这里
还有Freescience科研交流群
长按二维码加小秘书为好友
未经许可 不得转载长按二维码关注